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Abstract 

The  assumption  of  real-number  arithmetic,  which  is  at  the  basis  of  conventional 
geometric  algorithms,  has  been  seriously  challenged  in  recent  years,  since  digital  com¬ 
puters  do  not  exhibit  such  capability.  A  geometric  predicate  usually  consists  of  eval¬ 
uating  the  sign  of  some  algebraic  expression.  In  most  cases,  rounded  computations 
yield  a  reliable  result,  but  sometimes  rounded  arithmetic  introduces  errors  which  may 
invalidate  the  algorithms.  The  rounded  arithmetic  may  produce  an  incorrect  result 
only  if  the  exact  absolute  value  of  the  algebraic  expression  is  smaller  than  some  (small) 
e,  which  represent  s  the  largest  error  that  may  arise  in  the  evaluation  of  the  expression. 

The  threshold  e  depends  on  the  structure  of  the  expression  and  on  the  adopted  com¬ 
puter  arithmetic,  assuming  that  the  input  operands  are  error-free.  A  pair  (arithmetic 
engine, threshold)  is  an  arithinef.ir  filler.  In  this  paper  we  develop  a  general  technique 
for  assessing  the  efficacy  of  an  arithmetic  filter.  The  analysis  consists  of  evaluating 
both  the  threshold  and  the  probability  of  failure  of  the  filter.  To  exemplify  the  ap¬ 
proach,  under  the  assumption  that  the  input  points  be  chosen  randomly  in  a  unit  ball 
or  unit  cube  with  uniform  density,  we  analyze  the  two  important  predicates  ’’which- 
side"  and  ’insphere’\  We  show  that  the  probability  that  the  absolute  values  of  the 
corresponding  determinants  be  no  larger  than  some  positive  value  F,  with  emphasis 
on  small  V^  is  0(V)  for  the  which-side  predicate,  while  for  the  insphere  predicate  it  is 
Q(V^)  in  dimension  1,  O(V’i)  in  dimension  2,  and  0[V^  In  f )  in  higher  dimensions. 
Constants  are  small,  and  are  given  in  the  paper. 

1  Introduction 

The  original  model  of  Computational  Geometry  rests  on  real-number  arithmetic,  and  un¬ 
der  this  assumption  the  issue  of  precision  is  irrelevant.  However,  the  reality  that  computer 

•  This  work  was  partially  supported  by  ESPRIT  LTR  21957  (COAL)  and  by  the  U.S.  Army  Research  Ofi 
fice  under  grant  DAAH04‘96-t-0013.  This  work  was  done  while  0.  Devillers  was  visiting  Brown  University. 
This  document  is  being  simultaneously  released  by  l.N.R.I.A  as  Report  2971. 

^INRIA,  BP93,  06902  Sophia  Antipolis,  (France).  Dlivier.Devillersesophia.inria.fr 

^Brown  Univ.,  Dep.  of  Computer  Science,  Providence,  RI  02912-1910  (USA).  franco0cs.broHn.edu. 
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calculations  have  finite  precision  has  raised  an  increasing  awareness  of  its  effect  on  the  qual¬ 
ity  and  even  the  validity  of  geometric  algorithms  conceived  within  the  original  model,  in 
the  sense  that  algorithm  correctness  does  not  automatically  translate  into  program  cor¬ 
rectness.  In  recent  years  this  issue  has  been  amply  debated  in  the  literature  (see,  e.g., 
[BKM'*'95,  FV93,  Yap96]).  In  particular,  it  has  been  observed  that  while  some  degree  of 
approximation  may  be  tolerated  in  geometric  constructions,  the  evaluation  of  predicates 
(  the  ’’tests  ’’  carried  out  in  the  execution  of  programs,  -  such  as  which-side,  incircle, 
insphere  -)  must  be  exact  to  ensure  the  structural  (topological)  correctness  of  the  results 
[BMS94.  Yap96.  LPT96]. 

In  principle,  error-free  predicate  evaluation  is  achievable  for  error-free  input  operands, 
if  the  latter  are  treated  as  integers  and  the  arithmetic  is  carried  out  with  whatever  operand 
length  is  required  to  express  the  intermediate  results.  Such  safe  approach,  however, 
if  adopted  in  its  crudest  form,  would  involve  enormous  overheads  and  would  be  nearly 
impracticable,  since  the  execution  time  of  some  operation  (such  as  multiplication)  may 
increase  quadratically  with  the  length  of  the  representation. 

As  a  time-saving  alternative  to  exact  arithmetic,  it  has  been  customary  to  resort  to 
rounded  (approximate)  arithmetic  (for  example,  floating-point  arithmetic).  Such  practice 
can  be  modeled  as  follows.  Evaluation  of  a  predicate  P  typically  involves  computing  the 
value  /i  of  some  expression,  built  using  rational  operations.  The  value  of  /r  can  be  mapped 
to  one  of  three  values;  positive,  zero,  negative  (  referred  to  here  as  the  "sign”  of  /r), 
which  defines  the  predicate  P.  Let  S  denote  an  evaluator  for  P ,  and  let  denote  the 
numerical  value  computed  by  C.  In  general,  S  is  an  approximate  evaluator  of  /r,  so  that  its 
use  involve^  tlie  adoption  of  a  device,  called  a  certifier,  intended  to  validate  the  correctness 
of  the  evaluation.  The  pair  (evaluator,  certifier)  is  what  has  been  refered  to  as  a  filter  m 
the  literature  [F\'93.  MN94].  Typically,  the  certifier  for  f  compares  1^(5)1  with  a  fixed 
threshold  s(£)  >  0.  If  \n{£)\  >  e(S),  then  the  sign  of  n{£)  is  reliable.  Otherwise,  the 
certifier  is  unable  to  validate  the  result  and  we  have  a  failure  of  the  filter.  In  such  event 
recour'^e  to  a  more  powerful  filter  is  in  order.  This  suggests  the  need  to  develop  a  family 
of  filters  of  increasing  precision  (and  complexity),  to  be  used  in  sequence  until  failure  no 
longer  occurs.  The  last  item  of  this  sequence  is  the  exact  evaluator,  for  which  the  certifier 
is  vacuous  (i.e.,  s(^)  —  9).  Such  approach,  tvith  an  obvious  trade  off  bet^reen  efficacy  and 
efficiency,  embodies  the  notion  of  adaptive  precision. 

i.From  a  practical  standpoint  it  is  therefore  very  important  to  gauge  the  efficacy  of 
very  simple  filters,  that  is,  their  probability  of  success.  If  it  turn  out  that  if  a  filter  has 
a  high  probability  of  success,  then  recourse  to  a  more  time-consuming  filter  (or  exact 
computations)  will  be  a  rather  rare  event.  Of  course,  any  such  estimate  of  efficacy  rests 
on  some  arbitrary  hypotheses  on  the  o  priori  probability  of  problem  instances.  This  is 
an  important  caveat,  however,  under  reasonable  hypotheses  (uniform  distributions),  we 
submit  that  the  obtained  estimates  will  be  a  significant  contribution  to  the  assessment  of 
the  validity  of  such  approaches. 
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In  this  paper  we  analyze  two  specific  predicates,  “which-side”  and  “insphere”.  Both 
predicates  consist  of  computing  the  signs  of  appropriate  Sx6  determinants,  whose  entries 
are  specified  with  a  fixed  number  of  bits.  Depending  upon  the  adopted  evaluation  scheme 
(the  choice  of  an  equivalent  expression  for  a  given  function)  and  upon  the  precision  of 
the  input  operands  (for  example,  only  a  fixeddength  prefix  of  their  representation  may  be 
used  in  the  evaluation),  only  a  prefix  of  the  computed  value  is  reliable.  This  means  that 
if  the  absolute  value  of  the  determinant  is  above  a  known  threshold  e,  then  its  sign  is  also 
reliable. 

Our  objective  is  therefore  two- fold: 

1.  To  compute  the  value  of  the  threshold  e  for  a  given  determinantal  evaluation  tech¬ 
nique; 

2.  To  compute  the  probability  that  the  absolute  value  of  the  result  of  the  evaluation 
does  not  exceed  e,  i.e.,  the  probability  of  filter  failure. 

\Ne  recall  that  when  evaluating  the  determinant  of  a  (5  X  (5  matrix  A  (a  ^-determinant, 
for  short)  we  are  computing  the  signed  measure  of  a  hyperparallelepiped  defined  by  the 
8  vectors  corresponding  to  the  rows  of  A.  Since  each  of  the  components  of  these  vectors 
is  an  integer  in  the  range  -  1],  a  generic  vector  is  (applied  to  the  origin  and) 

defined  by  its  free  terminus  at  grid  points  in  a  (5-dimensional  cube  Cs  of  sidelength  2^ 
centered  at  the  origin. 

Our  probability  model  assumes  that  all  grid  points  within  Cs  have  identical  probability. 
Our  analysis  aims  at  estimating  (a  majorization  of)  the  distribution  of  the  volume  of 
the  hyperparallelogram  described  above.  To  obtain  the  desired  result,  we  introduce  some 
simplifications  consistent  with  the  objective  to  majorize  the  probability.  Specifically,  while 
our  assumption  is  a  discrete  uniform  distribution  within  we  begin  by  considering  a 
continuous  uniform  distribution  within  the  ball  Bs^,  the  <5-dimensional  ball  of  radius  1. 
The  obtained  results  are  then  used  to  bound  from  above  the  distribution  of  the  volume 
for  uniform  density  within  the  cube  Cs  (the  (5-fold  cartesian  product  of  interval  [—1, 1]), 
and  are  finally  extended  to  the  target  case  of  uniform  discrete  distribution.  We  shall 
recognize  that  the  initial  simplification  (uniform  density  in  Cs)  closely  approximates  the 
more  realistic  situation. 

The  paper  is  organized  as  follows.  We  begin  with  the  ”which-side”  predicate  (referred 
to  as  ’’determinant”),  i.e.,  in  Section  2  we  carry  out  the  probabilistic  analysis  in  Bs^  for 
S  =  1,2,3  and  arbitrary  6  (detailed  considerations  of  the  low-dimensional  cases  has  an 
obvious  pedagogical  motivation).  In  Section  3,  we  extend  the  results  from  the  continuous 
ball  to  the  discrete  cube.  In  Section  4  we  carry  out  an  analogous  analysis  for  the  ’’insphere” 
predicate,  which  illustrates  the  adverse  effect  of  dependencies  among  the  determinant 
entries.  Finally,  in  Section  7  we  evaluate  the  precision  of  determinant  evaluation  by 
recursive  expansion,  and  illustrate  the  efficacy /efficiency  tradeoff. 
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2  Probabilistic  analysis  within  unit  ball 

Throughout  this  section  we  adopt  the  following  notation:  We  let  Xi,i2)  •  • be  the 
coordinates  of  (J-space,and  pi,p2i  •  •  -Pi  be  <5  points  in  the  unit  ball  of  dimension  6.  We  also 
denote  by  \p\ ,  P2  •  •  -Pjl  the  absolute  value  of  the  determinant  defined  by  points  Pi,P2  •  ■  -Pj- 
This  quantity,  which  is  the  volume  of  the  hyperparallelogram  defined  by  the  origin  and  by 
points  P],P2  •  •  -Pii  s^lso  be  denoted  aj. 

We  begin  by  examining  in  some  detail  the  cases  of  low-dimension  determinants.  Since 
the  analysis  is  done  in  a  visualizable  geometric  setting  (^  <  3),  it  is  preparatory  to  the 
more  abstract  higher-dimensional  cases. 

2.1  1-and  2-deterininant 

Obviously,  if  pi  is  uniformly  distributed  between  -1  and  1,  then 

Prob{\pi\  <  R)  =  R 

Less  trivial  is  the  analy.sis  of  the  two-dimensional  case.  We  will  study  the  probability 
for  Ipi.pil  to  be  smaller  than  a  constant  A  when  pi  and  p2  are  distributed  uniformly  in 
the  unit  disk. 

Once  pi  is  chosen  (due  to  the  circular  symmetry,  pi  is  represented  by  a  single  parameter 
a-) ,  its  distance  from  the  origin),  P2  will  yield  an  area  between  (I2  and  (X2"b^®2  if  it  belongs 
to  one  of  the  two  strips  of  width  ^  depicted  in  Figure  1. 

We  then  have 


Fro&(|Pi,P2l  <  ^)  =  /  Froh(|pi,p2|  <  A  I  oi)p(ai)dai. 

Jo 

Since  the  density  function  of  oi  is  p(ai)  =  ^  =  2ai ,  and  the  density  of  02  conditional 
on  ai  is  p(a2  \  ai)da2  =  ^  ((^)  is  the  density  in  the  unit  disk,  and 

4^1  _  (szY  is  the  total  length  of  the  two  strips)  we  have 


Fro6(lp],p2|  <  A) 


min(>4,ai ) 


2(lidcLi  dcL2 


/*!  rminl 

Ja\—0  J 02^0 

Jo  (f 

I 


—da2 
0  TT 


02  arccos  ■ 


Ol 
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Figure  1:  For  the  analysis  of  the  2-determinant 


=  J  —da2  (^yjl  -  al  —  0,2  arccosa2^ 


8  3 

-  7^2 

TT  4 


arcsin  02 


-02  arccosa2 


0  2  ^ 

— A\/l  -  H —  arcsin  A -\ — arccos  A 

TT  TT  TT 


One  can  easily  verify  that  the  value  of  the  last  expression  is  1  for  A  =  1,  which  is  the 
maximum  attainable  value  for  the  area  of  the  parallelogram. 

2.2  3-determinant 

Again  we  assume  that  points  are  uniformly  distributed  in  the  unit  ball,  and  compute  the 
probability  that  this  volume  is  smaller  than  a  constant  V  >0. 

To  compute  the  volume  of  the  parallelepiped  defined  by  0,pi,p25  and  ps,  we  begin  by 
considering  the  parallelogram  defined  by  0,pi,  and  p2i  evaluate  its  area,  and  then  consider 
the  distance  of  pa  from  the  plane  containing  the  parallelogram. 

Distance  Opi  is  between  oi  and  oi  +  dai  if  pi  belongs  to  a  spherical  crown  of  thickness 
dai  and  area  47raj.  Therefore  the  probability  density  of  ai  is  p(ai)  =  ^47rai  =  SaJ  (note 
that  ^  is  the  density  in  the  unit  sphere). 

Once  Pi  has  been  chosen,  the  area  of  the  parallelogram  defined  by  0,pi  and  p2  is 
between  a2  and  02  +  da2  if  P2  belongs  to  the  crown  of  thickness  ^  of  a  cylinder  of  radius 
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Figure  2:  For  the  analysis  of  the  3-determinant 


^  whose  axis  contains  Opi  (see  Figure  2a).  Therefore  the  distribution  of  02  conditional 
on  ai  is  given  by 


p{a2\ai)da2 


(a) 


Finally,  once  and  p2  have  been  chosen,  the  volume  of  the  parallelepiped  is  between 
G3  and  0.3  +  daz  if  ps  belongs  to  one  of  the  two  spherical  slices  of  width  parallel  to  the 
plane  containing  0,pi,  and  p2,  and  at  distance  from  it  (see  Figure  2b).  Therefore  the 
distribution  of  03  conditional  on  03  is  given  by 


p(g3  I  02)^03  “  ^  1  ^ 


\02/  J  0-2  2a2  y  \0‘2/  / 


On  the  basis  of  this  analysis,  we  can  say 


rmin(V,02) 
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J d 02  =  0  03  =  0 
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f  —das  f  doi  ( ^dai  +  n— -  y/ai^  -  aa^—  -  ai  arcsin  —  -  2^  arcsin 

Jo  4  ^  \2  ai  ''  01  oi  01  «»/ 

jT  ^daa  lno3  -  ^oa  v/l  -  as*  +  (Soj^  +  i)arccosa3  -  2az^  j  ^arcsin^dai^ 
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.27  27  3 
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The  latter  integral  is  not  elementarily  computable.  Since  the  integrand  is  always  posi¬ 
tive,  so  is  the  integral.  To  neglect  it  corresponds  to  majorizing  the  probability.  Therefore 
we  write 


'27  TT 

Profc(|p],p2,p3|<  V)  <  -—y 
lb 


—  V  arcsin  V 
8 


—V^  -  —  InV-  —V®  arcsin  V 
8  8  4  4 


If  we  set  V'  =  1,  then  the  neglected  term  can  be  evaluated  exactly  to  ^  —  1  by 
exchanging  the  order  of  integration,  and  we  correctly  obtain  the  value  1  for  the  probability. 

The  preceding  analysis,  in  its  simplicity,  reveals  the  essential  items  for  the  evaluation 
of  the  relevant  conditional  probability  densities.  Specifically,  referring  concretely  to  the 
ca.se  (5  =  3,  due  to  the  assumption  of  uniform  distribution  of  the  points  in  the  unit 
sphere,  the  conditional  probability  density  p(oi|a,_i)da,-  of  a,-  given  a,_i,  i  =  2,3,  is 
proportional  (through  the  value  of  the  density  in  the  unit  sphere)  to  the  volume  of  some 
three-dimensional  domain.  The  latter  is  a  thin  crown  (of  thickness  of  a  three- 

dimensional  surface  53,,  which  is  the  locus  of  the  points  at  a  distance  between  and 

spanned  by  variables  xi,. .  .,a:,_i.  Surface  53,,-  has  a  very 
simple  structure.  Let  [u,  v)  be  a  pair  of  points  realizing  the  distance  with  u  G  ^i—i 
and  V  £  S3,,.  Point  v  belongs  to  the  boundary  of  a  (4  -  t)-dimensional  ball  of  radius 
of  which  point  u  is  the  center:  therefore  this  entire  boundary  belongs  to  53,,  (  in  the 
discussion  above,  this  boundary  consists  of  a  circle  for  i  =  2  and  of  two  points  for  i  =  3 
).  Moreover,  since  u  belongs  to  a  flat  any  translate  of  v  within  the  unit  sphere  in 

a  flat  parallel  to  also  belongs  to  53, j.  These  translates  form  the  intersection  of  the 
unit  sphere  with  a  flat  at  distance  from  the  center  of  the  sphere,  and  therefore  are  an 

_  l)-dimensional  ball  of  radius  ^1  -  (3^)  (in  the  discussionn  above  this  ball  consists 

of  a  segment  for  i  =  2  and  of  a  disk  for  i  =  3).  We  conclude  that  Ss^i  is  the  cartesian 
product  of  the  boundary  of  a  (4  -  i)-dimensional  ball  of  raxlius  (  the  "boundary” 

term)  and  of  an  {i  —  l)-dimensional  ball  of  radius  ^1  —  (57^7)  (ft®  "domain”  term). 
The  expression  fot  the  conditional  probability  density  consists  of  four  factors:  the  density 
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within  the  unit  sphere,  the  measure  of  the  boundary  term,  the  measure  of  the  domain 
term,  and  the  thickness  of  the  crown.  These  four  items  (in  the  order  given)  are  evidenced 
in  (a)  and  (b). 


2.3  Higher-dimensional  determinant 

We  now  extend  the  preceding  analysis  to  arbitrary  dimension  S.  If  we  assume  for  a  one¬ 
dimensional  volume  (a  distance)  the  conventional  degree  of  1,  then  the  volume  and  the 
surface  of  a  j-dimensional  domain  have  respective  degrees  j  and  j  —  1.  Let  Uj(r)  and 
Sj(r)  respectively  denote  volume  and  surface  of  ay-dimensional  ball  of  radius  r.  We  recall 
that[Ber87.  9.12.4.6] 


7r2 


u,(r)  =  —r*  for  i  even  u,(r)  = 


1-] 

2 


r'  for  i  odd. 


The  probability  density  p(ai)doi  =  Prob{ai  <  |Opi|  <  ai+dai)  is  obviously  given  by 
Referring  next  to  the  observations  at  the  end  of  the  preceding  subsection,  the 
conditional  probability  density  p(a,  |  a,_i)dai  of  a,  after  p\,P2i  •  •  -tPi-i  h^ve  been  chosen 
(conventionally  in  the  flat  described  by  coordinates  a:i,  X2,  ■ . JC.-i)  to  realize  the  value 
a,_] ,  has  the  following  expression: 


p(ai|a,_])dai 


i te)  ■“  (iRa' 


Therefore,  since  t’,(r)  =  and  5,(r)  =  i  •  Vi{l)  •  r*  we  have 
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p(ni )  JJp(a,  |a,-i)da: 


i=2 


In  the  rightmost  term  above  the  product  of  the  powers  of  the  a,’s  simplifies  to  so  that 
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The  last  expression  contains  a  constant  depending  only  on  <5,  which  will  be  denoted 


ks  =  - 


Ui(l) 


S-1 


Now  we  can  write  the  probability  for  the  absolute  value  of  the  determinant  to  be  no 
larger  than  V 


Prob{\p}  ^P2  ‘  ^  -  psl  <y)  —  ks 


r\  rai  /‘aa-2  /‘min(V,o^, 

t/ai  =  CK/a2=0  _  j  =(V  a£  =  0 


We  now  observe  that  the  integral 


/•min(V,ai_i) 

Jas=0 


das 


is  trivially  bounded  by  V .  Therefore  we  obtain  the  following  upper  bound: 
Prob(\p,.p,  ...nl<  V)  <k,vl^...p"  n  -  (^)  ) 


dot  j 


Since  the  latter  integral  does  not  depends  on  V,  its  value  is  another  constant,  which 
we  denote  Is  and  which  depends  only  on  S.  However,  when  V  =  1,  Equation  (c)  yields 
ProbHpj.p2  . .  .psi  <  1)  =  ksis+i  and,  since  Prob(lpi,P2  ■  ■  -Psl  <  1)  =  l,wegetli  = 
Finally,  we  conclude 


Prob(lpi,p2  ■  ■ -psl  <  V)  <  ksIsV 

vs(lp  * 

For  small  S,  the  bounds  are  given  below.  For  S  =  1,2,3  the  bounds  coincide  with  the 
values  previously  found,  that  is, 

CTi  <  V 
2^ 

(72  <  — « 2.5 

TT 
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3  7!-  r 

< 

2^ 

< 

213 

(T4 

3%2  ~  ■ 

< 

zW 

217  ' 

< 

224 

(^6 

5%3  ~  ' 

3  From  continuous  ball  to  discrete  cube 


3.1  From  continuous  ball  to  continuous  cube 

The  above  calculations  have  been  carried  out  for  points  uniformly  distributed  inside  the 
^-dimensional  ball  of  radius  1,  referred  to  djs>  Bs>  This  assumption  may  not  seem  to  model 
the  real  situation  for  two  reasons:  (i)  points  manipulated  by  computers  have  discrete 
rather  than  continuous  coordinates,  and  (ii)  points  are  more  reasonably  assumed  to  be 
uniformly  distributed  in  a  cube  than  in  a  ball  (as  in  the  case  when  each  coordinate  is 
independently  and  uniformly  selected). 

We  will  first  show  that  the  previous  result  relative  to  the  ball  Bs  induces  a  similar 
result  for  uniform  density  in  the  unit  cube  Cs  =  [—1, 1]^- 

Note  that  Cs  is  contained  within  a  5-dimensional  ball  of  radius  \/5,  denoted  y/SBs^ 
First,  we  consider  points  of  Cs  as  points  of  y/SBs  and  apply  a  homothety  with  a  factor 
thereby  obtaining 

Prob{\pi,p2. .  .ps\  <  ^  I  P.  €  y/6Bs)  =  Prob{\puP2  ■  ■ -Psl  <  V  \  pi  €  Bs)  < 


Next,  we  wish  to  restrict  the  points  to  belong  to  Cs,  i.e.,  we  consider  the  event  p,  € 
'/SBs-  i  =  as  the  union  of  the  event^p,  E  Cs,i  =  and  its  negation.  The 

probability  of  this  event  is  clearly  ( 


Prob{\puP2  •  •  -Psl  <  I  Pt  6  y/SBs)  >  ^  (i)v^)  ‘ ' 

which  gives  us  an  upper  bound  to  the  probablity  in  question.  Specifically 


Prob{\pi .  P2  .  •  -Pil  <  V’  I  Pi  G  Ci)  < 


( C^S  y  / 

vf  - 
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Practically,  for  small  values  of  S  we  get: 


V’l  = 
'tP2  = 

V’3  = 

t/>4  = 

= 

i>6  = 


1 


It  w  3.2 


27 

128' 

32 

8l’ 


-Tt 


-It 


21 


380 


9765625  12 

:7r 


402653184 
19683  15 


125000 


-Tt 


4.5 


23000 

•  10® 


The  previous  computation  can  probably  be  generalized  to  other  kinds  of  domains 
provided  that  the  ratio  between  the  volumes  of  the  inscribed  ball  and  of  the  circumscribing 
ball  is  bounded. 


3«2  From  continuous  cube  to  discrete  cube 

In  this  section  we  shall  discuss  why  the  obtained  results  for  continuous  density  of  points 
are  still  useful  for  discrete  probability,  i.e.,  when  points  belong  to  a  regular  grid  of  ^ 
points  inside  C$. 

Notice  that  the  preceding  results  are  clearly  incorrect  for  discrete  probabilities,  since 
in  the  latter  instance  the  probability  for  V  to  be  0  is  non-zero.  For  example,  in  two 
dimensions,  when  pi  is  chosen,  the  probability  that  p2  coincides  either  with  pi  or  with 
the  origin  is  2if,  and,  for  a  less  trivial  case,  the  probability  that  pi  has  both  coordinates 
smaller  than  5  is  5  and  thus  the  probability  that  p2  is  twice  as  large  as  is  However 
we  can  still  prove  that  Prob{\pi . .  .pi]  <  V)  =  OiV)  when  77  is  smaller  than  V . 

If  Pi  . .  .pi  is  a  set  of  points  in  Cs  whose  determinant  is  not  larger  than  F,  we  will  map 
it  to  a  the  nearest  set  of  grid  points  p'l . .  .p^  whose  determinant  is  not  too  large.  More 
precisely,  if  "P  =  (pi,  •  •  •  ,P^),  P'  =  (Pi,  •••iP's)  and  V  =  T'  -T  =  (dj, . .  .,<1^),  we  have 

|P'1  =  |P-HP1= 

where  7  is  a  nonempty  subset  of  {1,2..., 5}  and  |PP/|  is  the  determinant  obtained  by 
replacing,  for  each  i  6  7,  Pi  with  d,  in  [Pj.  The  above  result  follows  from  the  multilinearity 
of  the  determinant.  If  the  cardinality  of  7  is  j,  we  can  bound  from  above  the  absolute 
value  of  |PP;|  by  the  product  of  the  norms  of  its  vector.  Since  ||p,l|  <  y/6  and  l|dj||  <  y/S'^ 
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we  have  \V'Di\  <  By  grouping  the  (^)  terms  with  identical  value  of  j  we  get: 

\P'\-\V\  <y/f 

Referring  now  to  the  (J'^-dimensional  space  whose  points  are  the  sets  Pi...p5,  the 
determinant  \'P'\  is  no  larger  than  V  only  if  all  the  points  V  in  the  hypervoxel  (in  8^ 

dimensions)  have  determinant  \V\  no  larger  than  V  +  Since,  clearly,  a  random 

point  V  for  the  continuous  distribution  can  be  in  any  voxel  with  the  same  probability,  we 
conclude  that 

4  The  insphere  test 

In  the  preceding  analysis  of  the  ”which-side”  predicate,  the  points  defining  the  hyperpar¬ 
allelogram  were  assumed  to  be  independent  and  equally  distributed.  In  this  section  we 
consider  a  cause  for  which  there  exist  dependencies  among  the  coordinates  of  the  points: 
the  ’"insphere"  predicate.  This  predicate,  referred  to  as  <5-insphere  for  short,  tests  whether 
in  8  dimensions  the  origin  lies  inside  the  hypersphere  defined  by  {8  +  1)  other  arbitrary 
points  p,  =  z  =  1,  2, . . . ,  ^  +  1. 

It  is  well  known  that  the  (5-dimensional  insphere  test  is  embodied  in  the  sign  of  the 
determinant 


2^11 

X12 

Ijj  +  Xj2  +  .  . 

■  +  ^1S 

1 

2:21 

122 

.  X2J  "I"  X22  "1“  •  • 

■  +  ^ls 

1 

2^t5+2,l 

xs+2,2  ■ 

•  •  ^S+2,1  +  ^1+2,2  +  ‘ 

•'  +  ^S+2,S 

1 

Without  loss  of  generality,  one  of  these  points  (P(j^-2)  can  be  chosen  as  the  origin,  so 
that  the  above  determinant  becomes 


Xii 

X12 

•  ®11  *12  •  • 

■  +  ^u 

Ai  = 

X21 

X22 

x|i  +  X22  +  •  • 

■+4s 

2:^+1, 2  • 

••  ^S+hl+^Uh2  +  - 

•  •  +  *1+1, i 

These  8+1  points  are  assumed  to  be  evenly  distributed  in  the  unit  cube  C^. 
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4.1 


1-insphere  test 

We  can  model  the  problem  as  follows.  The  origin  O  and  u  define  the  1-dimensional  sphere; 
V  is  the  query  point.  Parameters  u  and  v  are  independent  and  uniformly  distributed  in 
[—1.1],  and  we  wish  to  evaluate  the  probability  of  the  following  event: 


l^il  <  ^ 


with  Aj  = 


u  V? 
V 


We  can  view  the  above  determinant  as  defined  by  two  points  («,  v)  and  (w^,  u^)  in  the 
u.  V  plane.  Clearlj'  the  choice  of  point  (u,  r)  completely  determines  the  determinant  value. 
Point  [u.  v)  is  uniformly  distributed  in  the  square  [—1, 1]  X  [—1, 1].  Since  uv{u  —  u)  is  the 
determinant  value,  the  determinant  is  null  on  the  three  lines  «  =  0,  v  =  0  and  u  =  v;  its 
values  are  symmetric  wdth  respect  to  the  line  u  =  -v  and  antisymmetric  with  respect  to 
u  =  V.  Therefore  for  any  value  of  0  <  A  <  2  it  is  sufficient  to  evaluate  the  probability  in 
the  quadrant  -u  <  v  <  u,0  <  u  <  1  (fully  shaded  quadrant  in  Figure  3)  and  multiply 
it  by  4.  In  the  upper  semiquadrant  (where  the  determinant  is  negative)  the  contour  lines 
have  equations: 


The  two  curves  for  fixed  A  join  with  a  common  vertical  tangent  at  the  point  {/j,,  |)  where 
=  (4.4)3  (notice  that  such  curves  exists  only  for  A  <  |).  In  the  lower  semiquadrant 
(where  the  determinant  is  positive)  contour  lines  have  equations: 


V  = 


u 

2 


V  =  —u 


and  intersect  the  line  v  =  -u  at  point  {u,  -v),  where  u  =  \{iA)'i  =  The  probability 
of  the  event  described  above  is  given  by  the  heavily  shaded  area  in  the  Figure  3  (in  fax:t  this 
area  should  be  multiplied  by  4  since  there  are  four  quadrants,  and  normalized,  dividing 
by  4,  since  4  is  the  area  of  the  square) .  The  area  is  also  1  minus  the  area  of  the  lightly 
shaded  region.  The  latter  area,  which  we  want  to  bound  from  below  is  given  by: 
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and  (for  A  <  ^) 


/i-44>1-44 


U' 


3  - 


SO  that 


J  >  C  udu- A  t  \du  +  I  udu-4A  [  \du 

Jif  Ji/ 

1  j/2  A  I  .  .  .  4A 

2  2^  i^  2  2  /i 

If  we  now  use  the  fact  the  fact  that  2u  =  /i  =  (4A)3 ,  then  we  obtain: 


1 7  ^/2  2  2 

Pro&(|Ai|  <  A)  <  -^^A3  -  5A  <  5.355^3 

A  direct  numerical  calculation  gives  Fro6(|Ai|  <  ^)  =  0.7,  while  the  value  of  the  above 
bounding  expression  for  the  same  value  of  A  is  0.85  (  an  excellent  agreement  considering 
the  rather  high  value  of  A).  For  A  >  |,  the  lightly-shaded  region  in  the  upper  semi¬ 
quadrant  disappears  yielding 

Prob{\A-i  \  <  A)  <  Sv^aI  -  4A  <  3.78At 

but,  considering  the  high  range  of  A,  this  expression  has  little  practical  interest. 

Alternatively,  we  may  consider  the  following  formulation.  Rather  than  lifting  point 
to  the  parabola  y  =  v'^u-  u^v  (for  fixed  positive  u) ,  we  lift  it  to  the  parabola  y  =  v'^-uv 
F{v)  (which  intersects  the  u-axis  at  0  and  u)  so  that 

Pro6(|A|  <  A)  =  Prob{\F{v)\  <  — ) 

'Vb 

This  formulation  will  be  useful  when  considering  the  multidimensional  case  in  Section 
4.2.2 


4.2  Higher  dimensional  insphere  test 

By  elementary  columns  operation  the  determinant  Aj  can  be  transformed  into  one  where 
the  last  column  has  zero  entries  except  in  the  last  row,  i.e.. 


X\1 

112 

Xis 

0 

X2l 

122 

X2S 

0 

XS\ 

XS2 

xss 

0 

XS+1,2  •  •  • 

w 
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Figure  3:  The  (u,  u)-region  defining  (|A|  <  A)  in  the  1-insphere  test 


The  determinant  of  the  intersection  of  the  first  6  rows  and  columns  of  the  above  matrix 
gives  the  signed  volume  v  of  the  hyperparallelogram  defined  by  the  first  S  points  and  W 
has  the  dimension  of  a  length  and  is  the  value  of  the  function 

,  ,,xs)  =  xj  +  xl  +  •  •  •  +  xj  -  cixi  -  . .  .csxs  (d) 

evaluated  at  point  In  the  expression  (d),  which  is  0  on  the  circumsphere  S  of 

pi, ....  PS.  is  the  center  of  the  above  circumsphere  and  =  4(^1  +  •  •  •  +  <^|) 

the  square  of  the  sphere  radius.  Notice  also  that  xs^i  =  F{xi^ . .  .,x^)  is  the  equation  of  a 
hyperparaboloid  H  in  {5  +  l)-dimensional  space,  so  that  . . . ,  =  F{ps.^i) 

is  the  signed  height  of  the  point  obtained  by  lifting  ps^i  from  the  hyperplane  =  0 
(to  which  it  belongs)  to  H.  Notice  that  hyperplane  =  0  contains  the  hypersphere  5. 

We  now  wish  to  bound  from  above  the  probability  of  the  event  |F(p5+i)|  <  V,  for 
some  constant  V, 

Assuming  as  usual  constant  density,  this  probability  is  the  volume  of  a  sphere  of 
radius  u"  =  y/V  +  x^  minus  the  volume  of  a  concentric  sphere  5'  of  radius  =  y/—V  + 
whenever  the  latter  is  defined.  (These  spheres  are  intersected  with  Cs  and  normalized  by  its 
volume.)  In  general  dimension,  since  the  above  radii  depend  upon  the  points  Pi,P2}  •  •  • 
the  evaluations  of  the  volumes  of  the  above  spheres  is  problematic.  However,  we  shall  show 
that  an  interesting  simplification  occurs  for  S  =  2. 

Below  we  shall  use  the  following  bounding  technique.  Let  As  be  expressed  as  the 


product  of  two  continuous  random  variables 

=  ab 

Given  a  constant  a  we  have: 

V 

Prob{ab  <V)  <  Prob{a  <  a  or  6  <  — -) 

V  y 

<  Probia  <  a)  +  Prohih  <  — )  -  Prob{a  <  a  and  b  <  —) 
—  \  ^  a  o; 


that  is: 


Prob{ab  <  V)  <  Prob{a  <  cv)  +  Prob{b  <  — ) 


(e) 


4.2.1  2-insphere  test 

In  this  case  the  volumes  of  S"  and  S'  are  respectively  ‘k[V  -f  v?)  and  tt[-V  +  “^))  ^ 
that,  if  S'  is  defined,  then  their  difference  becomes  7r2F,  otherwise  (i.e.,  when  v?  <  V) 
<  7r2V".  After  normalization  (since  4  is  the  measure  of  the  unit  square)  we  have 

Try 

Prob{\F{pz)\<V)<—. 


Given  that  IA2I  =  \P\P2\F{P?,)  inequality  (e)  becomes 

V 

Prob(\A2\  <V)  <  Prob{\piP2\  <  a)  -f  Pro6(|F(p3)|  <  — ) 

Try 

<  V>2a  +  — . 

From  Section  3  we  know  that  ^2  =  ^-  Choosing  a  suitable  value  of  a  we  have 

Prob{\A2\  <  y)  <  5.55v/y 


4.2.2  5-insphere  test  (5  >  2) 

As  for  the  two-dimensional  case,  we  shall  prove  that,  in  general,  the  probability  that 
|F(pi+i)|  is  no  larger  than  W  is  sufficiently  small. 

Again,  inequality  (e)  becomes: 

W 

Prob{\As\  <W)  <  Prob{\pi . .  .ps\  <  -^)  -f  Fro6(|F(p4+i)|  <  a)  (f) 
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We  know,  from  the  results  of  Section  3,  that 

W  W 
Prob{\pi . .  .ps\  <  — )  <  tps  — 

Therefore,  there  remains  to  bound  from  above  Pro6(|F(pi+i)|  <  o).  Recall  that 
|F(pf+i)l  is  the  distance  from  the  plane  Xi+i  =  0  of  the  point  lifted  to  the  hy¬ 
perparaboloid  U,  of  equation  F{p)  =  0,  which  intersects  xs+i  =  0  in  a  (5-dimensional 
sphere  passing  by  the  origin  with  with  center  q  =  (^, .  •  We  define  point  such 

that  it  lies  on  a  line  I  passing  by  O  and  q  and  such  that  length (p|^j, 9)  =  length (p^+i, 9) 
and  among  the  two  possible  choices,  we  select  the  one  closest  to  O  (see  Figure  4).  It  is 
immediate  that  |F(p<5+i)|  =  lF(pJ^j)|  and  that  length (OpJ_,.j)  <  \/S. 

Thus,  our  problem  is  reduced  to  a  one-dimensional  instance,  closely  related  to  the  one 
we  studied  in  Section  4.1  (here  u  =  21ength(09)  and  v  is  the  signed  value  of  length(OpJ+i)). 
There  are  however,  some  significant  differences.  There,  we  were  evaluating  the  probability 
of  the  event  |uF(u)|  <  V  ,  and  variables  n  and  v  were  uniformly  distributed  in  [-1,1]- 
Here,  on  the  other  hand,  we  wish  to  evaluate  the  probability  of  the  event  lF(u)|  <  o;,  ra¬ 
dius"  u  varies  between  0  and  oo,  and  “distance”  v  varies  between  -y/S  and  y/S.  Moreover, 
the  densities  pi  (u)  and  P2(i’)  of  u  and  u  respectively  are  not  constant;  however,  as  we  shall 
see.  they  are  bounded  by  constants  91  and  92- 

Next,  we  observe  that 


Prob{\F{ps+-i)\  <  o) 


/  /  pi(u)p2{v)Prob{\F(ps+i)\  <  01  \  u,v)dvdu 

Ju=0  Jv=-VS 

poo  p\/S 

2/  /  Piiu)p2{v)Prob{\F{ps+i)\<a\u,v)dvdu 

Ju=0  Jv=0 


Since  the  function  F{ps+\)  has  the  expression  -  uu,  the  right-hand-side  of  the  above 
equation  is  just  the  integral  of  Pi{u)P2iv)  on  the  domain(  shown  in  Figure  5)  bounded 
by  the  curves  —  uv  =  ±a.  This  integration  domain  can  be  split  into  four  subdomains 
A.B.C.D,  as  show-n  in  Figure  5.  This  split  depends  on  a  parameter  13  <  y/a  which  will  be 
chosen  later.  In  detail  we  have; 

•  Domain  A  is  defined  by  ^3  <  u  <  v/(^)  andu-f<u<t;-l-^.  Pi{u)p2{v)  <  9i92 
and  thus 

Jj  Pi{u)p2{v)dvdu  <  qiq2  ■— du  <  9192  ^In  <5 -4- 21n  — ^  a 

•  Domain  B  is  a  rectangle  defined  by  0  <  u  <  /?  and  0  j  <  u  <  (3+  pi{u)p2{v)  < 
QiQ]  and  thus 

rr  2a 

JJ^Pi{u)p2iv)dvdu  <  9i92/?-j^  <  29192a 
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•  Domain  C  is  defined  hy  0  <  v  <  /3  and  f3+^<u<v+^.  p2{v)  <  92  and  since 
/ p\(u)du  <  1  on  any  domain,  we  get 

JJ  p\{u)p2{v)dvdu  <  q2^ 

•  Domain  D  is  analogous  to  domain  C, 

JJ  Pi{u)p2{v)dvdu  <  92/3 

These  results  yield: 

Prob{\F{ps+i)\  <a)  =  29192®  ^In  (i  +  2 In  ^  +  2^  +  492/?.  (g) 
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To  obtain  values  for  q\  and  q2,  we  observe: 

1.  With  reference  topi  (ti),  consider  the  (^l)-dimensional  sphere  passing  by  0,p\. .  .ps-i 
in  the  flat  11  defined  by  these  points  (see  Figure  6  for  an  illustration),  and  let  t  be 
its  center.  Consider  the  family  T  of  ^-dimensional  spheres  passing  through  the  ori¬ 
gin  and  whose  centers  project  to  t  in  H.  Point  ps  will  determine  a  5-dimensional 
sphere  of  radius  between  u  and  u  -f  du  if  it  lies  between  a  pair  of  spheres  of  (or 
its  symmetric  pair  with  respect  to  11)  of  radii  u  and  u  -f  du.  The  volume  p\{u)du 
of  this  region  (normalized  by  the  volume  2^  of  Cs)  is  <  -^Asdu,  where  As  is  the 
maximum  of  the  measure  of  the  surface  of  a  5-dimensional  sphere  passing  by  the 
origin  and  intersected  with  Cs-  A  trivial  upper  bound  on  As  is  the  measure  252“^"^ 
of  the  surface  of  Cs-  Therefore 

Pi(«)<|.252^-'  =  25  4gi 

2.  With  reference  to  P2iv),  we  observe  that  once  ps  is  fixed,  the  choice  of  ps+i  will 
produce  a  value  between  v  and  v  +  dvii  ps+i  belongs  to  the  shaded  region  in  Figure 
7.  The  volume  of  this  region  is  dv  times  the  surface  of  the  sphere  of  radius  v  inside 
Cs-  This  surface  is  clearly  less  than  the  surface  of  Cs,  so  that  we  get 
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<f>3  «  70  X3  =  -100 
<t>4  «  408  X4  =  350 
(f>5  «  3970  X5  =  18000 
4>e  «  68500  X6  =  640000 

4.3  Discrete  distribution 

As  for  the  case  of  the  which-side  test,  also  for  the  insphere  test  we  can  map  each  point  to 
the  nearest  grid  point,  and  evaluate  the  ensueing  effect  on  the  determinant  to  be  computed. 
Generalizing  the  notation  of  Section  3.2,  we  have  that  the  norm  ||p'||  of  the  lifted  point  is 

bounded  by  and  that  ||(i|||  <  +  2^  <  +  Therefore,  grouping 

the  errors  term,  we  get 


and  considering  7]  small 


5  The  efficacy  of  arithmetic  filters 

To  complete  the  analysis,  in  this  section  we  wish  to  assess,  under  the  given  statistical 
assumptions,  the  probability  of  failure  of  a  given  arithmetic  filter  for  determinant  sign 
evaluation  (i.e..  the  probability  that  the  filter  is  unable  to  certify  the  correctness  of  the 
computed  sign).  This  probability  is  a  quantitative  measure  of  the  efficacy  of  the  filter.  If 
|pi  .  ..ps\  >  £6-  then  the  result  of  the  computation  is  reliable,  and  so  is  its  sign.  Plugging 
es  in  place  of  V  in  the  expression  above,  we  obtain  the  condition: 

P™6(|p....k|<.,|  O)  <  + 

where  we  have  used  the  result  (Section  3)  that  Prob{\pi . .  .ps\  <  V)  <  ‘tpsV .  The  param¬ 
eter  PS  introduced  here  is  therefore  the  sought  measure  of  filter  efficacy. 
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To  exemplify  this  approach,  we  shall  compute  ps  for  the  evaluation  of  determinants, 
for  the  case  where  the  coordinates  of  the  points  are  floating-point  numbers  in  the  interval 
[“1,1],  the  computations  are  carried  out  using  floating-point  arithmetic  with  6  bits  of 
mantissa,  and  the  determinant  is  evaluated  by  standard  expansion  with  respect  to  one  of 
its  columns  (recursive  evaluation). 

To  this  end,  it  is  necessary  to  compute  the  parameter  es-  We  introduce  the  following 
notation: 

•  5 [M.  77?]  denotes  the  set  of  numbers  whose  absolute  value  is  bounded  by  M  and 
whose  error  is  bounded  by  m.  Original  entries  belong  to  £[1,0]. 

•  M  denotes  2 

With  this  notation,  if  x\  E  £[Mi,r77i]  and  X2  E  £[M2,7?^2],  then  we  have: 

+  X2  €  S[My  d-  M2, 2”^"^  •  Ml  +  M2  +  mi  +  m2] 

'  X2E  S[M\  •  M2, 2~^“^  -  Ml  •  M2  +  mi  •  M2  +  m2  •  Ml] 

These  rules  express  the  mechanics  of  6-bit  mantissa  normalizing  floating-point  operations 
with  round-off.  After  transforming  variable  y  to  an  £[M,  m]  pair,  we  shall  express  the 
above  rules  as  an  arithmetics  on  such  pairs,  as  follows: 

1.  S\Mi.  777]]  -f"  £[M2»  77^2]  =  ^[Mi  +  M2, 2“^“*^  •  Ml  -f  M2  +  mi  +  m2] 

2.  £[M] .  77?i]  •  £[M2.  7772]  =  £[Mi  •  M2, •  Mi  •  M2  +  mi  •  M2  +  m2  •  Mi] 


The  results  given  below  are  obtained  in  an  Appendix  to  this  paper  by  the  mechanical 
application  (with  a  few  noted  exceptions)  of  the  above  two  rules  to  the  recursive  evaluation 
of  a  determinant.  Using  a  6-bit  mantissa,  the  results  are: 


•  €2  =  2-2-^ 

•  £5  =  576  •  2"*’ 

•  £8  =  226624  -  2"’’ 

•  £3  =  13  •  2-^ 

•  £6  =  3672  •  2"® 

•  £4  =  76  •  2"'' 

•  £7  =  27304  •  2-*’ 

For  example,  using  the  IEEE  norm  on  6  =  53  bits. 

we  can  therefore  estimate  the  cor- 

responding  probablility  of  failure.  The  pertinent  values  of  ,  Sy/o  and  ps  are  displayed 

below: 

•  £2  =  2.2-10-^® 

•  =  1.6-10-’® 

•  P2  =  1.2  - 10"’® 

•  £3  =  1.4- 10"’® 

•  =  8.7-10"’® 

•  P3  =  4.8  -  lO"’'* 
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•  £4  =  8.4  •  10'^® 


•  £5  =  5.7  • 

•  £6  =  8.3  •  10-^3 


=  7,1 .  io-’5 

<Jv^§  =  1.1 -10-12 


•  P4  =  5.9  •  10-12 

•  ps  =  3.0  •  10-^ 

•  pe  =  8.7  •  10-® 


Finally,  it  is  important  to  evaluate  the  efficiency  of  the  examined  filter.  We  observe 
here  that  a  careful  implementation  of  the  ’’recursive  evaluation  technique”  in  reality  should 
not  be  purely  recursive  (since  this  would  require  5!  operations).  In  fact,  the  recursive 
evaluation  involves  (^)  minors  of  dimension  i.  In  turn  each  such  minor  involves  (2i  —  1) 
arithmetic  operations  {i  multiplication  and  z  —  1  additions),  whose  operands  are  either 
minors  of  smaller  dimension  or  original  coefficients  (for  i  =  2,  where  the  recursion  stops). 
Thus,  the  total  number  of  operations  is: 


"<5  =  E  (  •)  -  i)  -  1)  =  (<^  -  1)(2'  -  1)  ~  <^2^ 

For  (5  <  7  we  obtain  rj  =  0,  r2  =  3,r3  =  14,  r4  =  45,  fs  =  124,  rg  =  315, r7  =  762  and 
rg  =  1785. 

It  must  be  pointed  out  that,  for  the  same  function  considered  above  (determinant 
value),  one  may  choose  alternative  evaluation  schemes  (corresponding  to  different  expres¬ 
sions  of  the  given  function,  such  as  other  expansion  rules,  Gaussian  elimination,  etc.) 
and/or  different  arithmetic  engines.  Each  such  choice  would  embody  a  filter,  whose  effi¬ 
cacy  and  efficiency  can  be  assessed  with  the  outlined  method. 


6  Summary  of  results  and  conclusions 

In  this  paper  we  have  developed  a  general  approach  to  the  cissessment  of  the  efficauzy  of 
arithmetic  filters,  under  some  reasonable  probability  assumptions.  As  an  important  exam¬ 
ple,  we  have  considered  the  efficacy  of  filters  for  the  evaluation  of  signs  of  determinants, 
both  for  a  case  where  all  entries  are  independent  (the  ”which-side”  predicate)  and  for 
a  case  where  dependencies  exist  (the  ’’insphere”  predicate).  This  analysis,  in  general, 
consists  of  two  parts. 

The  first  part  aims  at  computing  the  threshold  for  certification  by  the  filter  (i.e.,  the 
maximum  error  which  can  be  generated  by  the  evaluation  process),  and  does  not  rest  on 
any  assumption  about  the  distribution  of  the  data.  As  an  example,  we  have  carried  out 
this  threshold  analysis  for  the  so-called  recursive  evaluation  procedure,  which  computes  a 
determinant  by  expanding  it  with  respect  to  one  of  its  columns. 

The  second  part,  which  is  considerably  more  subtle,  aims  at  establishing  the  probability 
of  failure  of  the  filter,  i.e.,  the  probability  that  the  result  of  the  computation  falls  below 
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the  threshold.  This  analysis  rests  on  a  priori  assumptions  on  the  distribution  of  the  input 
data,  which  we  have  taken  as  uniform  within  their  representation  range,  and  has  been 
carried  out  for  the  two  important  geometric  tests  mentioned  above.  With  the  notations 
introduced  in  the  preceding  sections,  the  results  are  summarized  below. 


Prob(\p\,p2  •  ^  ^ps\  <  y  \  Pi  continuous  in  Bs)  < 

Prob{\pi,p2  -  -Psl  <V  \  Pi  continuous  in  Cs)  <  ^sV 

Prob{\puP2  . .  -Ps]  <V  \  pi  £  rpgridinCs)  <  tps^{V asn) 

Prob{\A}  I  <  I  Pi  continuous  in  Ci)  <  4AW^ 

Prob{\As\  <  W  I  p,'  continuous  in  Cs)  <  <i>sVW  ]n  —  d-XsVw 

Prob{\Ai\  <  W  \  Pi  €  rpgnd  in  Cs)  <  <l>s^/W  +  /3s^y^  ^  (dsy/p 


Below  we  repeat  for  synoptic  convenience  the  definitions  of  the  relevant  constants  and 
a  tabulation  of  their  values  for  small  8, 


<t>S 

XS 

Ps 


rVS-l{l)^ 


252 
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The  values  for  small  S  are 


s 

tl’S 

(l>s 

XS 

OiS 

0s 

1 

1 

1 

0.5 

2 

2 

2.5 

3.2 

5.6 

2 

18 

3 

5.3 

21 

70 

-100 

7.8 

200 

4 

10 

380 

408 

350 

32 

2800 

~ 

19 

23000 

3970 

18000 

140 

47000 

6 

35 

4.5  •  10^^ 

68500 

640000 

648 

900000 

The  above  values  are  tight  only  for  the  constant  a.  Due  to  data  interdependencies,  the 
analysis  of  the  insphere  predicate  is  considerably  more  involved  than  that  of  the  which-side 
predicate,  and  the  adverse  effect  of  the  dependencies  is  manifest  in  the  larger  values  of  the 
probability  of  failure. 

Such  analysis  is  particularly  valuable  for  estimating  the  time  required  to  test  the 
determinant  sign.  Under  the  given  probability  assumptions,  we  may  conclude  that  for 
small  dimension  (<  6)  straightforward  floating-point  filters  (i.e.,  floating-point  evaluators) 
are  extraordinarily  effective. 
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7  Appendix.  Error  evaluation 

D,  represents  a  generic  determinant  of  dimension  i,  and  z  a  generic  original  coordinate. 
It  must  be  pointed  out  that  the  recursive  evaluation  technique  yields  an  upper  bound  of 
(J!  for  Ds  because  dependencies  among  data  are  not  exploited.  However, smaller  values  of 
such  bound  are  known:  Since  the  determinant  value  is  bounded  by  the  product  of  the 
norms  of  its  S  components,  y/f  is  an  upper  bound  on  Ds,  which  is  attained  when  J  is  a 
power  of  two  (Hadamard  matrices).  For  small  values  of  5  the  following  bounds  have  been 
obtained  (  some  by  exhaustive  calculation):  Dj  <  1,  D2  <  2,  D3  <  4,  D4  <  16,  £>5  <  48, 
£>6  <  160.  Dt  <  576  and  Dg  <  4096.  Whenever  applicable,  we  shall  use  these  results 
below  to  obtain  tighter  estimates.  These  estimates  have  the  form 

Ds  6  €[Gs,Ss'\ 

where  zs  is  the  object  of  the  analysis  and  Gs,  the  largest  value  attainable  by  Ds,  is  only 
needed  to  carry  out  the  analysis. 

1.  D\  G  £[1.  0].  by  definition. 

2.  £>2  G  £[2.2“*’+’].  In  fact  £>2  =  z-£>i  +  z-£»i,  z  G  £[1,0]  and  £>i  G  £[1,0].  Therefore 

£[1.0] -£[1,0] +  £[1,0] -£[1,0]  =  £[l,2-'’-i]  +  £[l,2-'‘-J]  =  £[2,2-*’+i] 

3.  £>3  G  £[4, 13  •  2“**].  In  fact,  £>3  =  (z  •  £>2  +  a;  •  £>2)  +  ^  ‘  •C>2-  Therefore: 

£>3  G  (£[1, 0]  •  £[2, 2-'‘+i]  +  £[1, 0]  •  £[2, 2-''+!])  +  £[1, 0]  •  £[2, 2-*’+^] 

=  £[4,2-^-i.4  +  6-2-*>]  +  £[2,3-2-^] 

=  £[4, 2-*’-^  •  4  +  11  •  2-^]  =  £[4, 13  •  2"'’] 

where  we  have  used  the  fact  £>4  <  4. 

4.  £>4  G  £[16, 19  •  2"*’+^].  The  result  is  obtained  by  appljdng  the  previous  rules  and 
z  G  £[1, 0],  £>3  G  £[4, 13  •  2“*’]  to  the  following  evaluation  scheme: 

D4  =  (z  •  D3  +  Z  •  1)3)  +  (z  ■  D3  +  Z  •  D3) . 
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5.  I>5  e  ^[48, 129  •  2"^+2j  The  result  is  obtained  by  applying  the  previous  rules  and 
X  G  5[1, 0],  D4  G  ^^[16, 19  •  2“*’'^^]  to  the  following  evaluation  scheme: 

Dz  =  ({x-  D4  +  X  •  D4)  +  X  •  D4)  +  {x-  D4  +  X’  D4). 
making  also  use  of  the  fact  Ds  <  48. 

6.  De  G  ^^[160, 459  •  2“*^^].  The  result  is  obtained  by  applying  the  previous  rules  and 
X  E  S[l,  0],  r>5  G  5[48, 125  •  2"*’+^]  to  the  following  evaluation  scheme: 

De,  =  (((2;  •  Z>5  4- 1  ■  D5)  +  (x  •  D5  +  a;  ■  D5))  ■{■{x  ■  D^  +  x-  D5)) 
making  also  use  of  the  fact  De  <  160. 

7.  D7  G  5[576, 3413  •  2“*’+^].  The  result  is  obtained  by  applying  the  previous  rules  and 

X  G  ^^[1.  0],  De  G  £^[160,  459  •  to  the  following  evaluation  scheme: 

Dt  =  {{{x  •  Dq  +  X  ■  De)  +  [x  ■  De  +  X  •  De))  +  {{x  ■  De  +  x  ■  De)  +  x  ■  De)) 
making  also  use  of  the  fact  Dj  <  576. 

8.  De  G  £’[4096,3541-2“'’+®].  Again,  the  result  is  obtained  by  applying  the  previous 
rules  and  x  G  ^[l,  0],  D7  G  5[576, 3413  •  2“'’+®]  to  the  following  evaluation  scheme: 

Dg  —  ({{x  ■  D7  +  X  ■  D7)  +  (a:  •  D7  +  z  •  P?))  +  ((2;  •  P?  +  ^  •  P?)  +  (z  •  P7  +  z  •  P7))) 

making  also  use  of  the  fact  Pg  <  4096. 

REM.4RK.  If  the  original  values  do  not  belongs  to  £^[1,0]  but  to  £^[l,e]  and  e  is  small 
enough,  then  the  final  error  on  P,  must  be  increased  by  f!e.  c  small  enough  means  that 
the  error  on  M  does  not  affect  the  value  M,  in  particular  if  all  the  manipulated  sets 
are  of  the  form  S[M,  m]  with  M  =  M  +  m,  then  adopting  a  first-order  approximation  is 
legitimate. 
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